{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению\n",
    "<center>Автор материала: Каменев Михаил Сергеевич."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Калибровка вероятностей в ML</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Калибровка? Какая калибровка?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Приветствую тебя, дорогой читатель! Прости меня за этот неформальный стиль, надеюсь, он тебя не бесит :p Не хотелось бы сразу давать тучу теоретических выкладок, доказывающих, что некоторые алгоритмы в методе predict_proba возвращают вовсе не proba, а что-то непонятное. Поэтому, рассмотрим одну очень простую, но чрезвычайно полезную для понимания, задачу из теории связи. Давай мы с тобой разберемся в том, как работает демодулятор!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Если ты далёк от радиофизики, ничего страшного, я тоже от неё далёк :D. Но я уверен, что у тебя есть телефон, а значит ты каждый день демодулируешь. Когда базовая станция передает на твой телефон сигнал, он представляется в виде каких-то точек на комплексной плоскости. Например, следующий набор точек соответствует созвездию QPSK."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from scipy.stats import multivariate_normal\n",
    "from sklearn.metrics import roc_auc_score\n",
    "from sklearn.preprocessing import StandardScaler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qpsk_points = np.array([[-1,-1],[-1,1], [1,-1], [1,1]])\n",
    "sns.set_style(\"darkgrid\")\n",
    "sns.regplot(x=qpsk_points[:,0], y=qpsk_points[:,1],fit_reg=False);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Соответственно, базовая станция выбирает одну из этих точек и передает на телефон. Всё было бы хорошо, если бы не шум в канале связи. К сожалению, телефон получает не точку $x$, а $y = x + \\mathcal{N}(0, \\sigma ^ 2)$. Выглядит это как-то так:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(4444)\n",
    "points_count = 100\n",
    "\n",
    "transmitted_points = np.random.randint(2,size=(points_count, 2))\n",
    "classes = 2 * transmitted_points[:,0] + transmitted_points[:,1]\n",
    "transmitted_points = transmitted_points * 2 - 1\n",
    "noise = np.random.multivariate_normal([0,0], [[0.1,0],[0,0.1]], size=points_count)\n",
    "received_points = transmitted_points + noise\n",
    "\n",
    "data = pd.DataFrame(data=np.hstack((received_points, classes.reshape(-1,1))), columns=['x','y', 'target'])\n",
    "sns.lmplot(x='x', y='y', data=data, hue='target', fit_reg=False);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Поэтому демодулятор вычисляет вероятность того, что базовая станция передала одну из возможных точек при условии принятой точки $y$: $P(x = x' | y)$. Естественно, вычисляется это по формуле Байесса:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ P(x = x' | y) = \\frac{P(y | x = x')P(x')}{P(y)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Вероятности по всем возможным четырем значениям в дальнейшем будут нормироваться, чтобы их сумма была равна единице. При этом дробь $\\frac{P(x')}{P(y)}$ будет принимать одинаковые значения для всех четырех значений (предполагаем, что P(x') = 1/4 для всех $x'$, то есть все точки равновероятны). Поэтому, её можно не вычислять: она всё равно уйдет после нормировки (слава нормировке!). Из всей этой неясной писанины следует одно: достаточно посчитать $P(y | x = x')$ для всех возможных значений $x'$ и отнормировать эти значения так, чтобы они в сумме давали 1. Тем самым мы получим честные значения условных вероятностей!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сгенерим выборку из 100000 точек и посчитаем честные вероятности способом, описанным выше. Затем обучим логистическую регрессию, добавив полиномиальные признаки. Из четвертой лекции мы знаем, что без них логистическая регрессия в данной задаче будет работать очень плохо. Затем сравним честные вероятности с теми, что вернёт метод predict_proba класса LogisticRegression (Disclaimer: всё будет отлично!). И в конце сделаем всё тоже самое, но для случайного леса (Disclaimer: всё будет очень плохо)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Генерируем зашумленные точки\n",
    "np.random.seed(4444)\n",
    "points_count = 100000\n",
    "\n",
    "transmitted_points = np.random.randint(2,size=(points_count, 2))\n",
    "classes = 2 * transmitted_points[:,0] + transmitted_points[:,1]\n",
    "transmitted_points = transmitted_points * 2 - 1\n",
    "noise = np.random.multivariate_normal([0,0], [[1,0],[0,1]], size=points_count)\n",
    "received_points = transmitted_points + noise\n",
    "\n",
    "# Считаем условные вероятности\n",
    "probs = np.zeros((points_count, 4))\n",
    "all_possible_points = np.array([[-1, -1], [-1, 1], [1, -1], [1,1]])\n",
    "for row_idx, received_point  in enumerate(received_points):\n",
    "    for col_idx, point in enumerate(all_possible_points):\n",
    "        probs[row_idx, col_idx] = multivariate_normal.pdf(received_point, mean=point, cov=[[1,0],[0,1]])\n",
    "\n",
    "\n",
    "# Делаем нормировочку\n",
    "probs_norm = probs / np.repeat(np.sum(probs, axis=1).reshape(-1,1),4,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.model_selection import train_test_split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Добавляем полиномиальные признаки и делим нашу выборку на обучающую и тестовую\n",
    "X = np.hstack((received_points, received_points ** 2, (received_points[:,0] * received_points[:,1]).reshape(-1,1)))\n",
    "y = classes\n",
    "X_train, X_test, y_train, y_test, probs_train, probs_test = train_test_split(X,y,probs_norm, test_size=0.3,shuffle=True, random_state=4444)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Обучаем логистическую регрессию\n",
    "logit = LogisticRegression(random_state=4444, multi_class='multinomial', solver='lbfgs')\n",
    "logit.fit(X_train, y_train)\n",
    "probs_logit = logit.predict_proba(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# смотрим на ошибку\n",
    "np.mean(np.abs(probs_logit - probs_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видно, что в среднем логистическая регрессия ошибается на 0.002. На мой взгял, очень даже неплохой результат. Посмотрим теперь на случайный лес!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "random_forest = RandomForestClassifier(random_state=4444, n_estimators=100)\n",
    "random_forest.fit(X_train[:,[0,1]], y_train) # отбрасываем полиномиальные признаки\n",
    "probs_forest = random_forest.predict_proba(X_test[:,[0,1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(np.abs(probs_forest - probs_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "0.08. Как минимум, гораздо больше чем 0.002. На всякий случай, посмотрим на accuracy_score для обоих алгоритмов."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(accuracy_score(y_test, np.argmax(probs_logit, axis=1)))\n",
    "print(accuracy_score(y_test, np.argmax(probs_forest, axis=1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как вариант, можно попробовать случайный вес с деревьями максимальной глубины 2. Точность в этом случае станет одинаковой, но с вероятностями всё будет чуть хуже."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Дорогой читатель, тут мы подходим к главной причине написания этого тьюториала: случайный лес, svm, байессовские методы косячат с вероятностями. Сильно косячат. Поэтому было предложенно два метода (Platt's scaling и Isotonic Regression), чтобы как-то это дело подправить. И еще:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/andrew_ng_dts.jpg\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Некоторые вводные"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Начнем c того, что в документации sklearn есть отличная <a href='http://scikit-learn.org/stable/modules/calibration.html'>статья</a> на рассматриваемую тему. Для меня её было немножко тяжеловато читать и в ней совсем не рассказывается о том, как работать с методами калибровки. Однако там есть очень красивые картинки, среди которых присутствуют кривые калибровки.\n",
    "\n",
    "Предположим, что модель предсказывает вероятность принадлежности объекта к целевому классу в 0.5. Рассмотрим все объекты, которым модель предсказала вероятность в промежутке [0.49, 0.51]. Тогда среди рассматриваемых объектов, около половины должно принадлежать целевому классу (если алгоритм выдает честные вероятности, конечно).\n",
    "\n",
    "Для начала, оставим в нашем примере только две точки: [1,1] и [-1,-1]. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Генерируем зашумленные точки\n",
    "np.random.seed(4444)\n",
    "points_count = 100000\n",
    "\n",
    "transmitted_points = np.random.randint(2,size=(points_count, 1))\n",
    "classes = transmitted_points[:,0]\n",
    "transmitted_points = np.repeat(transmitted_points * 2 - 1,2,1)\n",
    "noise = np.random.multivariate_normal([0,0], [[1,0],[0,1]], size=points_count)\n",
    "received_points = transmitted_points + noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Метод для построения кривой калибровки\n",
    "def calibration_curve(y_pred, y_true):\n",
    "    sorted_idx = np.argsort(y_pred)\n",
    "    y_pred_sorted = y_pred[sorted_idx]\n",
    "    y_true_sorted = y_true[sorted_idx]\n",
    "    bins = np.linspace(0,1,11)\n",
    "    bin_values = np.zeros((bins.size - 1, 2))\n",
    "    for idx in range(1,bins.size):\n",
    "        from_prob = bins[idx - 1]\n",
    "        to_prob = bins[idx]\n",
    "        positive_count = 0\n",
    "        total_count = 0\n",
    "        for prob, y in zip(y_pred_sorted,y_true_sorted):\n",
    "            if prob >= from_prob:\n",
    "                total_count += 1.\n",
    "                if y == 1:\n",
    "                    positive_count += 1.\n",
    "            if prob > to_prob:\n",
    "                break\n",
    "        bin_values[idx - 1, :] = [(from_prob + to_prob) / 2, positive_count / total_count if total_count != 0 else 0]\n",
    "    return bin_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.hstack((received_points, received_points ** 2, (received_points[:,0] * received_points[:,1]).reshape(-1,1)))\n",
    "y = classes\n",
    "X_train, X_test, y_train, y_test = train_test_split(X,y, test_size=0.3,shuffle=True, random_state=4444)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logit = LogisticRegression(random_state=4444)\n",
    "logit.fit(X_train, y_train)\n",
    "probs_logit = logit.predict_proba(X_test)[:,1]\n",
    "random_forest = RandomForestClassifier(random_state=4444, n_estimators=100)\n",
    "random_forest.fit(X_train[:,[0,1]], y_train) # отбрасываем полиномиальные признаки\n",
    "probs_forest = random_forest.predict_proba(X_test[:,[0,1]])[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ideal_curve = np.linspace(0,1,21)\n",
    "fig, ax = plt.subplots(figsize=(12,12))\n",
    "#ax = plt.plot(x=ideal_curve, y=ideal_curve, ax=ax, linestyles='dashed')\n",
    "#ax = plt.plot(x=bins[:,0], y=bins[:,1], ax=ax, color='red')\n",
    "plt.plot(ideal_curve, ideal_curve, '--')\n",
    "\n",
    "logit_curve = calibration_curve(probs_logit, y_test)\n",
    "plt.plot(logit_curve[:,0], logit_curve[:,1], 'o-',color='r')\n",
    "forest_curve = calibration_curve(probs_forest, y_test)\n",
    "plt.plot(forest_curve[:,0], forest_curve[:,1], 'o-',color='g')\n",
    "plt.legend(['ideal', 'logit', 'random forest'])\n",
    "plt.xlabel('Probability')\n",
    "plt.ylabel('Fraction of positive')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Чудеса! Может показаться, что всё хорошо и где-то в начале была допущена грубая ошибка. Но давайте, сначала, рассмотрим реальные данные. Например, датасет с информацией по оттоку клиентов."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# код сворован из тетрадки с четвертой лекции\n",
    "\n",
    "data = pd.read_csv('../../data/telecom_churn.csv').drop('State', axis=1)\n",
    "data['International plan'] = data['International plan'].map({'Yes': 1, 'No': 0})\n",
    "data['Voice mail plan'] = data['Voice mail plan'].map({'Yes': 1, 'No': 0})\n",
    "\n",
    "y = data['Churn'].astype('int').values\n",
    "X = data.drop('Churn', axis=1).values\n",
    "X_train, X_test, y_train, y_test = train_test_split(X,y, test_size=0.3,shuffle=True, random_state=4444)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logit = LogisticRegression(random_state=4444)\n",
    "random_forest = RandomForestClassifier(random_state=4444, n_estimators=100)\n",
    "logit.fit(X_train,y_train)\n",
    "random_forest.fit(X_train, y_train)\n",
    "probs_logit = logit.predict_proba(X_test)[:,1]\n",
    "probs_forest = random_forest.predict_proba(X_test)[:,1]\n",
    "\n",
    "ideal_curve = np.linspace(0,1,21)\n",
    "fig, ax = plt.subplots(figsize=(12,12))\n",
    "plt.plot(ideal_curve, ideal_curve, '--')\n",
    "\n",
    "logit_curve = calibration_curve(probs_logit, y_test)\n",
    "plt.plot(logit_curve[:,0], logit_curve[:,1], 'o-',color='r')\n",
    "forest_curve = calibration_curve(probs_forest, y_test)\n",
    "plt.plot(forest_curve[:,0], forest_curve[:,1], 'o-',color='g')\n",
    "plt.legend(['ideal', 'logit', 'random forest'])\n",
    "plt.xlabel('Probability')\n",
    "plt.ylabel('Fraction of positive')\n",
    "plt.show()\n",
    "\n",
    "print('logit roc-auc score: %f' % roc_auc_score(y_test,probs_logit))\n",
    "print('random forest roc-auc score: %f' % roc_auc_score(y_test,probs_forest))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как известно, чудес не бывает =). Несмотря на то, что случайный лес по качеству оказался лучше логита, по вероятностям логит выглядит симпатичнее. Задача калибровки сводится к тому, чтобы приблизить хорошую модель, но с плохими вероятностями, к идеальной кривой. Тем самым мы не потеряем в качестве, а с другой стороны будем получать честные вероятности."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Platt's scaling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Данный алгоритм был предложен в 1999 товарищем по фамилии Platt в <a href='http://citeseer.ist.psu.edu/viewdoc/summary?doi=10.1.1.41.1639'> статье </a> Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods. Из названия видно, что алгоритм изначально был предложен для улучшения svm. Однако он работает для любой модели, в том числе и для случайного леса. \n",
    "\n",
    "Идея алгоритма проста до безобразия: \n",
    "- создаем отложенную выборку;\n",
    "- обучаем модель на обучающей выборке;\n",
    "- делаем предсказания на отложенной выборке;\n",
    "- обучяем логистическую регрессию на этих предсказаниях;\n",
    "- profit!\n",
    "\n",
    "В статье есть маленькое уточнение для целевой переменной в четвертом пункте. Обычно метки класса равны -1 и 1. Однако Platt рассматривает значения $\\frac{N_+ + 1}{N_+ + 2}$ и $\\frac{1}{N_- + 2}$. В статье вывод этих формул отсутствует и, по крайней мере мне, не совсем понятно, почему рассматриваются именно такие значения.\n",
    "\n",
    "Для наглядности, попробуем реализовать данный метод."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1) Создаем отложенную выборку\n",
    "X_train_platt, X_valid_platt, y_train_platt, y_valid_platt = train_test_split(X_train,y_train, test_size=0.3,shuffle=True, random_state=4444)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2) Обучаем случайный лес\n",
    "random_forest.fit(X_train_platt, y_train_platt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3) Делаем предсказание для отложенной выборки\n",
    "probs_forest_valid = random_forest.predict_proba(X_valid_platt)[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4) Обучаем логистичекую регрессию на предсказаниях леса\n",
    "# в идеале, тут бы применить scale из статьи, но тогда logit выбросит ошибку \n",
    "# так что, если строго, это не совсем platt scaling =)\n",
    "logit_platt = LogisticRegression(random_state=4444)\n",
    "logit_platt.fit(probs_forest_valid.reshape(-1,1), y_valid_platt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 5) profit!\n",
    "probs_forest_platt = logit_platt.predict_proba(random_forest.predict_proba(X_test)[:,1].reshape(-1,1))[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ideal_curve = np.linspace(0,1,21)\n",
    "fig, ax = plt.subplots(figsize=(12,12))\n",
    "plt.plot(ideal_curve, ideal_curve, '--')\n",
    "\n",
    "logit_curve = calibration_curve(probs_logit, y_test)\n",
    "plt.plot(logit_curve[:,0], logit_curve[:,1], 'o-',color='r')\n",
    "forest_curve = calibration_curve(probs_forest, y_test)\n",
    "plt.plot(forest_curve[:,0], forest_curve[:,1], 'o-',color='g')\n",
    "forest_curve_platt = calibration_curve(probs_forest_platt, y_test)\n",
    "plt.plot(forest_curve_platt[:,0], forest_curve_platt[:,1], 'o-',color='k')\n",
    "plt.legend(['ideal', 'logit', 'random forest', 'random forest (our platt)'])\n",
    "plt.xlabel('Probability')\n",
    "plt.ylabel('Fraction of positive')\n",
    "plt.show()\n",
    "\n",
    "print('logit roc-auc score: %f' % roc_auc_score(y_test,probs_logit))\n",
    "print('random forest with platt roc-auc score: %f' % roc_auc_score(y_test,probs_forest_platt))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видно, что после калибровки (черная кривая) вероятности стали чуть лучше в сравнеии с \"чистым\" случайным лесом (зеленая кривая)\n",
    "\n",
    "Известно, что изобретать велосипед не очень хорошо, поэтому воспользуемся реализованным в sklearn классом CalibratedClassifierCV. Его конструктор принимает на вход следующие параметры:\n",
    "\n",
    "- base_estimator - модель, которую будем калибровать\n",
    "- method - тут возможно два варианта: 'sigmoid' - это рассматриваемый platt scaling и 'isotonic' - это isotonic regression, к ней мы еще вернемся\n",
    "- cv - параметры кросс-валидации. По умолчанию данные бьются на три фолда. Можно указать количество фолдов с помощью целого числа. Можно передать объект-генератор. Если коротко, то это обычный cv, как и во многих других моделях. Есть одно но: можно передать в качестве cv строчку \"prefit\". Тогда вся выборка будет рассматриваться как отложенная (модель то уже обучена)\n",
    "\n",
    "Проделаем шаги 1-5, но уже с помощью CalibratedClassifierCV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.calibration import CalibratedClassifierCV\n",
    "\n",
    "X_train_platt, X_valid_platt, y_train_platt, y_valid_platt = train_test_split(X_train,y_train, test_size=0.3,shuffle=True, random_state=4444)\n",
    "random_forest.fit(X_train_platt, y_train_platt)\n",
    "calibrator = CalibratedClassifierCV(base_estimator=random_forest, method='sigmoid', cv='prefit')\n",
    "calibrator.fit(X_valid_platt,y_valid_platt)\n",
    "probs_forest_platt_calibrator = calibrator.predict_proba(X_test)[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ideal_curve = np.linspace(0,1,21)\n",
    "fig, ax = plt.subplots(figsize=(12,12))\n",
    "plt.plot(ideal_curve, ideal_curve, '--')\n",
    "\n",
    "logit_curve = calibration_curve(probs_logit, y_test)\n",
    "plt.plot(logit_curve[:,0], logit_curve[:,1], 'o-',color='r')\n",
    "forest_curve = calibration_curve(probs_forest, y_test)\n",
    "plt.plot(forest_curve[:,0], forest_curve[:,1], 'o-',color='g')\n",
    "forest_curve_platt = calibration_curve(probs_forest_platt, y_test)\n",
    "plt.plot(forest_curve_platt[:,0], forest_curve_platt[:,1], 'o-',color='k')\n",
    "forest_curve_platt_calibrator = calibration_curve(probs_forest_platt_calibrator, y_test)\n",
    "plt.plot(forest_curve_platt_calibrator[:,0], forest_curve_platt_calibrator[:,1], 'o-',color='y')\n",
    "plt.legend(['ideal', 'logit', 'random forest', 'random forest (our platt)', 'random forest (sigmoid)'])\n",
    "plt.xlabel('Probability')\n",
    "plt.ylabel('Fraction of positive')\n",
    "plt.show()\n",
    "\n",
    "print('logit roc-auc score: %f' % roc_auc_score(y_test,probs_logit))\n",
    "print('random forest with sigmoid callibration roc-auc score: %f' % roc_auc_score(y_test,probs_forest_platt_calibrator))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Тут мы должны извлечь еще один урок: не стоит реализовать то, что уже реализовано кем-то =). Но желтая кривая уже выглядит заметно лучше зеленой."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Isotonic regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Начнем сразу с примера"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_iso, X_valid_iso, y_train_iso, y_valid_iso = train_test_split(X_train,y_train, test_size=0.3,shuffle=True, random_state=4444)\n",
    "random_forest.fit(X_train_iso, y_train_iso)\n",
    "calibrator = CalibratedClassifierCV(base_estimator=random_forest, method='isotonic', cv='prefit')\n",
    "calibrator.fit(X_valid_iso,y_valid_iso)\n",
    "probs_forest_iso_calibrator = calibrator.predict_proba(X_test)[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ideal_curve = np.linspace(0,1,21)\n",
    "fig, ax = plt.subplots(figsize=(12,12))\n",
    "plt.plot(ideal_curve, ideal_curve, '--')\n",
    "\n",
    "logit_curve = calibration_curve(probs_logit, y_test)\n",
    "plt.plot(logit_curve[:,0], logit_curve[:,1], 'o-',color='r')\n",
    "forest_curve = calibration_curve(probs_forest, y_test)\n",
    "plt.plot(forest_curve[:,0], forest_curve[:,1], 'o-',color='g')\n",
    "forest_curve_platt = calibration_curve(probs_forest_platt, y_test)\n",
    "plt.plot(forest_curve_platt[:,0], forest_curve_platt[:,1], 'o-',color='k')\n",
    "forest_curve_platt_calibrator = calibration_curve(probs_forest_platt_calibrator, y_test)\n",
    "plt.plot(forest_curve_platt_calibrator[:,0], forest_curve_platt_calibrator[:,1], 'o-',color='y')\n",
    "forest_curve_iso_calibrator = calibration_curve(probs_forest_iso_calibrator, y_test)\n",
    "plt.plot(forest_curve_iso_calibrator[:,0], forest_curve_iso_calibrator[:,1], 'o-',color='c')\n",
    "plt.legend(['ideal', 'logit', 'random forest', 'random forest (our platt)', 'random forest (sigmoid)', 'random forest (isotonic)'])\n",
    "plt.xlabel('Probability')\n",
    "plt.ylabel('Fraction of positive')\n",
    "plt.show()\n",
    "\n",
    "print('logit roc-auc score: %f' % roc_auc_score(y_test,probs_logit))\n",
    "print('random forest with isotonic callibration roc-auc score: %f' % roc_auc_score(y_test,probs_forest_platt_calibrator))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Тут видна одна маленькая особенность изотонической регрессии: если данных мало, она переобучается. В документации приводится цифра <<1000, у нас же объектов было 700. Но всё таки, что же это такое \"изотоническая регрессия\"? На просторах интернета есть замечательная картинка:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/isotonic_example.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Изотоническая регрессия стремится аппроксимировать данные с помощью не одной линии, а нескольких, но при этом они должны образовывать неубывающую функцию. Очевидно, что данный тип регрессии, в общем случае, сможет подстроиться под данные лучше, чем обычная. Всё упирается в сложность. <a href=\"http://stat.wikia.com/wiki/Isotonic_regression#Pool_Adjacent_Violators_Algorithm\">Утверждается</a>, что в общем случае сложность обучения составляет порядка $n^4$, где $n$ - число объектов. Однако в нашем случае базовый алгоритм (случайный лес) хорошо ранжирует объекты, хоть и с неправильными вероятностями. Если $x_i, x_j$ - это вероятности для двух объектов с точки зрения случайного леса, а $y_i, y_j$ - реальные вероятности, можно смело утверждать, что если $x_i < x_j$, то и $y_i < y_j$. А если реальные вероятности отранжированны, есть простой <a href=\"http://stat.wikia.com/wiki/Isotonic_regression#Pool_Adjacent_Violators_Algorithm\">алгоритм</a>, который обучается за линейной время."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Дабы всё таки убедиться, что изотоническая регрессия работает, рассмотрим dataset из десятого домашнего задания: предсказание задержки авиарейса."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# код частично сворован из 10ого дз\n",
    "train = pd.read_csv('../../data/flight_delays_train.csv')\n",
    "X_train, y_train = train[['Distance', 'DepTime']].values, train['dep_delayed_15min'].map({'Y': 1, 'N': 0}).values\n",
    "\n",
    "X_train_part, X_test, y_train_part, y_test = train_test_split(X_train, y_train, test_size=0.3, random_state=17)\n",
    "\n",
    "scaler = StandardScaler()\n",
    "X_train_part = scaler.fit_transform(X_train_part)\n",
    "X_test = scaler.transform(X_test)\n",
    "# обучаем логистическую регрессию\n",
    "logit.fit(X_train_part, y_train_part)\n",
    "probs_logit = logit.predict_proba(X_test)[:, 1]\n",
    "# обучаем случайный лес\n",
    "random_forest.fit(X_train_part, y_train_part)\n",
    "probs_forest = random_forest.predict_proba(X_test)[:, 1]\n",
    "\n",
    "# делаем отложенную выборку для калибровки\n",
    "X_train_calibr, X_valid_calibr, y_train_calibr, y_valid_calibr = train_test_split(X_train_part,y_train_part, test_size=0.3,shuffle=True, random_state=4444)\n",
    "#platt scaling\n",
    "random_forest.fit(X_train_calibr, y_train_calibr)\n",
    "calibrator = CalibratedClassifierCV(base_estimator=random_forest, method='sigmoid', cv='prefit')\n",
    "calibrator.fit(X_valid_calibr,y_valid_calibr)\n",
    "probs_forest_platt_calibrator = calibrator.predict_proba(X_test)[:,1]\n",
    "# isotonic regression\n",
    "calibrator = CalibratedClassifierCV(base_estimator=random_forest, method='isotonic', cv='prefit')\n",
    "calibrator.fit(X_valid_calibr,y_valid_calibr)\n",
    "probs_forest_iso_calibrator = calibrator.predict_proba(X_test)[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ideal_curve = np.linspace(0,1,21)\n",
    "fig, ax = plt.subplots(figsize=(12,12))\n",
    "plt.plot(ideal_curve, ideal_curve, '--')\n",
    "\n",
    "logit_curve = calibration_curve(probs_logit, y_test)\n",
    "plt.plot(logit_curve[:,0], logit_curve[:,1], 'o-',color='r')\n",
    "forest_curve = calibration_curve(probs_forest, y_test)\n",
    "plt.plot(forest_curve[:,0], forest_curve[:,1], 'o-',color='g')\n",
    "forest_curve_platt_calibrator = calibration_curve(probs_forest_platt_calibrator, y_test)\n",
    "plt.plot(forest_curve_platt_calibrator[:,0], forest_curve_platt_calibrator[:,1], 'o-',color='y')\n",
    "forest_curve_iso_calibrator = calibration_curve(probs_forest_iso_calibrator, y_test)\n",
    "plt.plot(forest_curve_iso_calibrator[:,0], forest_curve_iso_calibrator[:,1], 'o-',color='c')\n",
    "plt.legend(['ideal', 'logit', 'random forest', 'random forest (sigmoid)', 'random forest (isotonic)'])\n",
    "plt.xlabel('Probability')\n",
    "plt.ylabel('Fraction of positive')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('logit roc-auc score: %f' % roc_auc_score(y_test,probs_logit))\n",
    "print('random forest roc-auc score: %f' % roc_auc_score(y_test,probs_forest))\n",
    "print('random forest with platt scaling roc-auc score: %f' % roc_auc_score(y_test,probs_forest_platt_calibrator))\n",
    "print('random forest with isotonic callibration roc-auc score: %f' % roc_auc_score(y_test,probs_forest_iso_calibrator))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Все рассмотренные модели имеют схожее качество, однако применение изотонической регрессии позволило лучше всего приблизить вероятности к идеальной прямой. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Заключение "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сегодня мы узнали (а может, дорогой читатель, ты уже это знал), что predict_proba не всегда возвращает честные вероятности. Однако существует два метода, которые позволяют приблизить ответы модели к честным вероятностям. Platt's scaling прост и работает в любой ситуации, в то время как изотоническая регрессия имеет свойство переобучаться на малом количестве данных. Но если данных достаточно, предпочтительным является использование именно изотонической регрессии. В данном тьюториале основной упор был сделан на практику, потому что есть множество статей, которые в деталях описывают, как эти методы работают:\n",
    "\n",
    "- википедия (<a href=\"https://en.wikipedia.org/wiki/Isotonic_regression\"> изотоническая регрессия </a> и <a href=\"https://en.wikipedia.org/wiki/Platt_scaling\">Platt's scaling </a>)\n",
    "- еще одна <a href=\"http://stat.wikia.com/wiki/Isotonic_regression\">вики</a> с подробным описанием алгоритма оптимизации изотонической регрессии\n",
    "- Статья Платта <a href='http://citeseer.ist.psu.edu/viewdoc/summary?doi=10.1.1.41.1639'>Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods</a>. Сам метод там описан хорошо\n",
    "-  <a href=\"http://scikit-learn.org/stable/modules/calibration.html\">Страница</a>, посвященная калибровке вероятностей на sklearn. Очень много хороших картинок, которые я тут уже дублировать не стал."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
